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Abstract 

Source confusion defines a practical depth to which to take large-area extragalactic surveys. 3D imaging 
spectrometers with positional as well as spectral information, however, potentially provide a means by 
which to use line emission to break the traditional confusion limit. In this paper we present the results 
of our investigation into the effectiveness of mid/far infrared, wide-area spectroscopic surveys in breaking 
the confusion limit. We use SAFARI, a FIR imaging Fourier Transform Spectrometer concept for the 
proposed JAXA-led SPICA mission, as a test case. We generate artificial skies representative of 100 
SAFARI footprints and use a fully-automated redshift determination method to retrieve redshifts for both 
spatially and spectrally confused sources for bright-end and burst mode galaxy evolution models. We find 
we are able to retrieve accurate redshifts for 38/54% of the brightest spectrally confused sources, with 
continuum fluxes as much as an order of magnitude below the 120 /zm photometric confusion limit. In 
addition we also recover accurate redshifts for 38/29% of the second brightest spectrally confused sources. 
Our results suggest that deep, spectral line surveys with SAFARI can break the traditional photometric 
confusion limit, and will also not only resolve, but provide redshifts for, a large number of previously 
inaccessible galaxies. To conclude we discuss some of the limitations of the technique, as well as further 
work. 
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1. Introduction 

The cosmic infrared background (CIB) peaks at ~ 150 
/an and comprises the total infrared (IR) emission from 
all sources in the sky (eg. Dole et al. 2001, Elbaz et 
al. 2002), integrated over all time. It has been found 
to contain as much energy as the combined optical/UV 
extragalactic background, suggesting that half of all light 
emitted by stars and active galactic nuclei (AGN) is ab- 
sorbed by dust before we are able to observe it in the 
optical (Hauser & Dwek 2001). Locally, the IR output 
of typical galaxies is only one third of their optical out- 
put (Soifer & Neugebauer 1991), which implies strong 
evolution in the IR properties of galaxies as one moves to 
high redshift. The CIB has been well-studied using many 
instruments including ISOCAM (ISO), MIPS (Spitzer), 
ISOPHOT (ISO) and SCUBA(JCMT) at 15, 24, 160 and 
850 /zm (Elbaz et al. 2002, Papovich et al. 2004, Juvela 
et al. 2000 and Smail et al. 2002 respectively). Although 
the peak of the CIB lies at ~ 150 /zm, it has yet to be re- 
solved into individual sources at these wavelengths: our 
understanding of the make-up of the CIB therefore relies 
on extrapolation. 

Observations in the mid and far-infrared (MIR and FIR; 
typically defined as lying in the wavebands 5-30 and 30- 
1000 /zm respectively) from ground based telescopes are 
difficult if not impossible because of the high opacity of the 
Earth's atmosphere at these wavelengths. Ground based 



telescopes are only sensitive to narrow wavebands in the 
MIR/FIR where the atmosphere's transmission is higher, 
thus wide band observations in the MIR/FIR must typ- 
ically be made from space. As a result, FIR telescopes 
are smaller in diameter than their optical counterparts 
(up to a few meters), with a resulting angular resolu- 
tion that is low compared to both optical and radio tele- 
scopes/interferometers. A direct consequence of this low 
angular resolution is that completely resolving the CIB 
into discrete sources in the FIR is very difficult, if not 
impossible, because of source confusion. 

Source confusion may be defined as the degradation of 
the quality of photometry of sources clustered on a scale 
to the order of the telescope beam size (eg. Scheuer 
1957). The confusion limit sets the useful depth to which 
large-area extra-galactic surveys should be taken. For ex- 
ample, the Herschel mission (Pilbratt 2004), successfully 
launched in May 2009, has a 3.5 m diameter mirror which 
realizes an angular resolution of 8" at 120 /zm. At these 
wavelengths, the confusion limit for such a mirror is esti- 
mated to be around ~5 mJy (eg. Dole et al. 2004, Jeong 
et al. 2006). Surveys at 24 /zm suggest that at these flux 
levels one will only be able to resolve at most ~ 50% of 
the CIB (Dole et al. 2004). It is possible to reduce the 
confusion limit through making observations with a larger 
diameter mirror. However, due to the practical limitations 
of high angular resolution FIR imaging there is a limit on 
how much we can reduce confusion noise. 
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One way to break the confusion limit makes use of the 
extra dimension of wavelength, to which one has access in 
spectroscopic surveys. Discrete sources can be identified 
by relatively bright, narrow-band emission lines: thereby 
allowing redshifts to be determined. A preliminary study 
to explore the efficacy of blind, wide area spectroscopic 
surveys in resolving FIR sources is described in Clements 
et al. (2007). In their work an artificial 'sky' was pop- 
ulated using template FIR spectral energy distributions 
(SEDs) of a selection of different types of galaxy, to which 
were added FIR emission lines of strengths derived from 
ISO-LWS observations (eg. Negishi et al. 2001). The 
sources were rcdshifted and assigned luminosities accord- 
ing to the evolutionary models of Pearson (2005), Pearson 
et al. (2007) and Pearson & Khan (2009). Observations 
of the 'sky' were made using the instrumental parameters 
(eg. sensitivity/noise levels, spectral resolution, field of 
view (FoV), beam size) of SAFARI, a FIR imaging Fourier 
Transform Spectrometer concept for the proposed JAXA- 
lcd SPICA (Space Infrared Telescope for Astronomy and 
Astrophysics) (Swinyard et al. 2008) mission. It will of- 
fer the large FoV and high spectral resolution required to 
break the confusion limit using spectroscopy. According 
to its current specifications SPICA will have a 3.5 m di- 
ameter mirror, and therefore will be subject to the same 
confusion noise as Hcrschcl. The primary mirror will, how- 
ever, be cooled to <6 K and so will offer a great leap in sen- 
sitivity over Herschel. SAFARI will cover the waveband 
35 to 210 /jm with varying resolution, including R^IOOO - 
(at 120 /im, A A = 0.176 /mi, when run in SAFARI'S higher 
resolution mode) - which matches the typical width of an 
cxtragalactic MIR/FIR emission line. 

Estimates of source redshift were made by hand by lo- 
cating the position of the strongest emission line in each 
spectrum. The strongest lines typically observed in the 
FIR are the [OI] and [CII] lines at 63.18 and 157.74 /mi 
respectively. If one assumes a typical dust temperature of 
35 K, then the strongest lines shortward and longward of 
the peak of the SED will be [OI] and [CII] respectively 
Beyond z = 2.5 these lines are shifted out of SAFARI'S 
observable waveband, therefore the redshifts of more dis- 
tant sources than this are irretrievable. If these two lines 
were the only ones present in the FIR then by compar- 
ing the evaluated source redshifts with the model input 
redshifts, one can assess the efficiency of this blind-line 
method. It was found that when looking at a patch of 
simulated 'sky' equal to one SAFARI FoV, it was possi- 
ble to retrieve accurate redshifts for sources with 120 /im 
continuum fluxes as much as a factor of -~10 below the 
traditional continuum confusion limit. Sources with 120 
/im flux S'i20/im > 1 mJy and at redshifts z < 2.5 were 
retrieved with 100% accuracy. 

The use of blind spectral line surveys to resolve FIR 
sources is not without its own limitations and type of con- 
fusion. Line, or spectral confusion occurs when multiple 
sources are observed in a single telescope beam: the spec- 
tra from two or more objects are effectively scrambled, and 
it can become difficult to determine which lines are emit- 
ted by which objects. As a result, source redshifts become 



hard to extract. The work described in Clements et al. 
(2007) made use of model spectra with FIR emission lines 
only. To assess the true viability of using spectral line sur- 
veys to break the confusion limit requires the inclusion of 
MIR emission lines in the 'sky' model, as sources will, in 
general, have both FIR and MIR emission lines. Inclusion 
of these shorter wavelength lines will enable the recov- 
ery of sources with redshifts of z > 2.5, beyond which the 
[OI] and [CII] emission lines at 63.18 and 157.74 /im, re- 
spectively, are shifted out of the SAFARI waveband. By 
including MIR lines however, one increases the problem 
of line confusion, and so assigning lines to individual, but 
spatially unresolved, sources becomes more problematic. 

In this paper we examine a much larger model 'sky' 
populated with more realistic template spectra with both 
FIR and MIR emission lines and employ a new automated 
method of evaluating source redshifts in a time efficient 
manner. We also implement a method to extract the red- 
shifts of multiple sources clustered in a single spatial bin. 
Through the implementation of this method we investi- 
gate how effectively we can break the traditional photo- 
metric confusion limit. In sections 2.1 and 2.2 we describe 
the generation of the artificial sky used to test the source 
recovery technique, which in turn is outlined in section 2.3. 
This is followed in section 3 by a quantitative assessment 
of the retrieval and error rates of the redshifts output, and 
a discussion of the results in section 4. 

2. Extracting Redshifts From an Artificial Sky 

Deep, blind-field imaging spectroscopy has the poten- 
tial to enable discrete sources at fluxes below the tradi- 
tional continuum confusion limit to be extracted. Using 
imaging spectroscopy, it is theoretically possible to extract 
redshifts for all line emitting sources present in the instru- 
ment's FoV, allowing all sources to be discretely resolved. 
In this work this is done using an automated redshift de- 
termination algorithm. In order to best test the efficiency 
(number of sources for which we determine redshifts, in- 
versely weighted by how many sources we inaccurately 
determine redshifts for) of this method we generate an ar- 
tificial sky in the form of a datacube, populated with re- 
alistic spectra taken from nearby galaxies, and redshifted 
according to the bright-end and burst mode galaxy evolu- 
tion models of Pearson (2005), Pearson et al. (2007) and 
Pearson & Khan (2009). Running the program through 
the datacube we compare the fitted redshifts with their in- 
put values in order to determine the method's precision, 
accuracy and efficiency 

2.1. Generating an Artificial Sky 

We create an initial set of two 1 square degree arti- 
ficial skies with 1" spatial pixels, populated by galax- 
ies with redshifts, spectral types and 40 /im continuum 
fluxes based on the bright-end and burst mode evolution- 
ary models of Pearson (2005), Pearson et al. (2007) and 
Pearson & Khan (2009). We do not include any physi- 
cally based spatial distribution modelling (eg. clustering, 
etc.) and in this work the skies are populated with sources 



3 



uniformly distributed in random positions. It should also 
be noted that we do not include any cirrus contribution 
and we assume SAFARI to be background limited. These 
act as our master 'skies' which we later crop and re-bin to 
create the SAFARI footprints. 

The bright-end and burst mode evolutionary models are 
backward evolution formulations where observed galaxy 
source counts are used to constrain the model parameters. 
The model components consist of a luminosity function 
to represent the number density of sources as a function 
of luminosity, a library of spectral energy distributions 
(SED) to model the extragalactic source population emis- 
sion as a function of wavelength and an assumption on the 
type- dependent evolution of the extragalactic population 
(in luminosity and number density). Both evolutionary 
models utilize the IRAS infrared local luminosity func- 
tion defined at 60/xm (Saunders et al. 2000) or 12/xm 
(Rush ct al. 1993) for the galaxy and AGN (Seyfert) 
populations respectively. Although various other, more 
recent luminosity functions are available, the IRAS func- 
tions has the advantage of being defined at or around the 
peak of the population emission spectrum and are free 
of contamination by mid-infrared features. The 60/im 
galaxy luminosity function is segregated into cool (nor- 
mal galaxy) and warm (star-forming) components, de- 
fined by IRAS colours where cool 100 /im/60/im cirrus-like 
colours (Efstathiou & Rowan-Robinson 2003) represent 
the normal quiescent galaxy population and the warmer 
100 /jm/60 /im colour component is representative of star- 
forming galaxies with activity increasing as a function of 
luminosity for M82-like Starburst Lju < lO n L , lumi- 
nous (LIRG) L IR > lO n L & ultraluminous (ULIRG) 
Lin > 1O 12 L0 infrared galaxies. Thus the model frame- 
work includes five general evolutionary population classes 
(Normal, Starburst, LIRG, ULIRG, AGN) of extragalac- 
tic object defined by luminosity, colour and subsequent 
evolution. 

The bright-end model assumes evolution in the galaxy 
population in both the density and luminosity for sources 
modeled by simple power laws of the form f(z) = 
(1 + z) k , where k is the type dependent evolutionary 
strength parameter. The evolutionary model is an up- 
dated framework of that first presented in Pearson & 
Rowan- Robinson (1996), with modest starburst galax- 
ies rather than ULIRGs dominating in 15 /tm selected 
source counts. The burst mode evolutionary model pre- 
dicts that the upturn of emission at 15 /mi (Elbaz et al. 
1999) and peak at 24 /xm (Papovich et al. 2004) is to 
due the emergence of a new population of L/ULIRGs. 
The original burst mode evolutionary model presented by 
Pearson (2001) has power law evolution similar to the 
bright-end model for the starburst and AGN sources and 
an initial violent exponential evolutionary phase of the 

form, f{z) = 1 + f.exp[— ^2^ ] ; from z = to z p = 1, 
where where k & a are the type dependent evolutionary 
strength parameters, followed by a power law evolution- 
ary phase for the L/UILRGs. Both the bright-end and 
burst mode evolutionary models have non-evolving nor- 



mal galaxy populations. Throughout this work, values re- 
garding the different evolution models are written in the 
form bright-cnd(burst mode). 

Each of the five general evolutionary population compo- 
nents (Normal, Starburst, LIRG, ULIRG, AGN) are rep- 
resented by a small set of galaxy spectral energy distribu- 
tions from the libraries of Efstathiou & Rowan-Robinson 
(2003), Efstathiou et al. (2000) and Efstathiou & Rowan- 
Robinson (1995) for the normal, starburst, L/ULIRG and 
AGN types respectively. The selected SEDs are represen- 
tative of the SED libraries from which they have been 
drawn and have been shown to be consistent with the 
colours of sources detected in the European Large Area 
ISO Survey (Oliver et al. 2000). Rowan- Robinson et 
al. (2004) showed that the ISO infrared galaxy pop- 
ulation could indeed be divided into four general spec- 
tral classes; normal quiescent galaxies; starburst (M82- 
likc) galaxies; luminous and ultraluminous infrared galax- 
ies and AGN. Rowan-Robinson ct al. (2004) found 
that the normal, quiescent population of ISO galaxies 
were well modelled with the templates of Efstathiou & 
Rowan- Robinson (2003) with far/mid-infrared ratios of 
vS u (100fj.m)/vS u (12(im)~6-7 whilst in contrast, IRAS 
galaxies were often modelled with templates with ratios 
of ^5 (Rowan-Robinson & Crawford 1989). Therefore, 
the normal galaxy component consist of 2 SEDs each of 
^S y (100/im)/VS , ! ,(12/im)'-^5.8 & 6 referred to as a normal 
and cold normal type respectively. For the star form- 
ing (starburst, L/ULIRG) population we have selected 
7 SEDs from the template libraries to ensure a varia- 
tion in the the mid-infrared features in an attempt to 
avoid artifacts caused by a particular choice of SED for 
all sources. The SEDs in the template libraries have a 
broad correlation between the model optical depth and 
the galaxy luminosity. 2/7 of these SEDs are selected 
for the starburst component (ry ~ 50 of which one is a 
model for the archetypal star- forming galaxy M82). 3/7 
of these SEDs are selected with increasing optical depths 
(ry ~ 50 — 100) for the LIRG component (assuming a cor- 
responding increase in luminosity for each SED of 10 n Lo, 
1O 115 L and a colder SED of L > lO n - 5 L referred to 
as a cold-LIRG) and the remaining 2/7 SEDs selected 
for the ULIRG component correspond to the best tem- 
plate model fits for the archetypal ULIRGs Arp220 (cold 
ULIRG) and Mk231 (hot ULIRG) respectively. Given the 
relatively featureless infrared spectra of AGN, the AGN 
component SED corresponds to a single tapered disc dust 
torus model. Emission lines arc then added to each of 
the template spectra, taken from ISO-LWS observations 
of nearby galaxies (eg. Negishi et al. 2001). 

MIR emission line strengths are taken, where possible, 
from the same sources as are used for the FIR emission 
lines. In some cases data were not currently available at 
these shorter wavelengths in which case sources with sim- 
ilar FIR characteristics as the original template spectra 
are used. The galaxies from which we take the MIR emis- 
sion line strengths are listed in table 1. The model spectra 
are shown in figure 1 and are added to the master 'sky' 
in spectral resolution of AA = 0.08 /im (R = 2500 at 200 
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Component Galaxy Type Line Template Reference 

Normal Cold/Normal NGC 7331 Smith et al. 2004 

Starburst M82/Starburst M82 Forster Schreiber el al. 2001 

LIRG lO n / n - 5 L /Cold LIRG NGC 253 Sturm el al. 2000 

ULIRG Hot/Cold ULIRG Arp 220 Sturm et al. 1996 

AGN Seyfert 1 Mrk 1014 Annus et al. 2004 

AGN Seyfert 2 NGC 1068 Lutz et al. 2000 



Table 1. Sources used for the addition of MIR lines to the template spectra that populate the data cube. 




Fig. 1. Template spectra which were used to populate the datacube. From the top downwards are the SED templates for; normal 
cold, normal, starburst M82, starburst, 10 n L Q LIRG, 10 11 ' 5 Lq LIRG, cold LIRG, hot ULIRG, cold ULIRG, Seyfert 1 and Seyfert 2. 



/im) in the waveband from 1 to 400 /im. 

The flux and redshift distributions for each SED type 
for the bright-end and burst mode evolution models are 
shown in figure 2. These figures show the distribution of 
SED type with 120 /im flux and redshift over a 1 square 
degree region of sky. In these regions of sky there are a to- 
tal of 25596(38975) sources with Sno^m. > 0.342 mJy (la, 
lOhr sensitivity of SAFARI). Of these sources 33(58)% 
have Si20nm > 1 rnJy. The burst mode evolution model in- 
cludes more high flux sources as well as more high redshift 
sources than the bright-end model. The burst mode model 
also includes the more extreme luminous/ultraluminous 
and Seyfert spectral types, which tend to dominate at 
higher rcdshifts. 100(75)% of sources with Si2o^ m > 1 mJy 
have redshifts z < 2.5, however setting Si20/mi > 0.342 mJy 
these values drop to 95(74)%. 

2.2. Generating Datacubes/ Source Catalogs 

SAFARI (Swinyard et al. 2008) will have a 2'x2' FoV 
with a diffraction-limited angular resolution of 8" at 120 
/im. The instrument will cover the waveband from 35 
to 210 /im at varying spectral resolution (R^IOOO at 120 
/im). However, in this work we match our data cubes 
to the original specifications of SAFARI, taking a FoV of 
128" xl28" and a waveband covering from 30 to 210 /im. 
Thus to generate datacubes representative of the sky as 



would be seen by SAFARI we take 128" x 128" sections 
of our artificial skies and re-bin to a spatial resolution 
of 8" (to simplify we have assumed a wavelength inde- 
pendent spatial and spectral resolution). This re-binning 
allows for the possibility of two or more sources in a sin- 
gle spatial bin. In these cases we refer to the brightest 
source with Svzo^m > 0.342 mJy as the primary source 
and the second brightest >5i20^m > 0.342 mJy as the sec- 
ondary source. Applying the angular resolution of SPICA 
means that 6(8)% of all pixels have two or more sources 
present. If there are multiple sources present in a spa- 
tial pixel only the two brightest will be investigated as 
other sources will be too faint relative to the primary and 
secondary sources to detect. We create our spectra by 
cropping the waveband of the artificial sky to between 30 
and 210 /im with and smoothing them to a fixed resolu- 
tion of AA = 0.176 /im. In this work we use the original 
wavelength range SAFARI rather than its current speci- 
fication, thus wc have a slightly larger waveband to pick 
up emission lines from. We then add Gaussian noise with 
a standard deviation of a = 0.342 mJy and a zero mean 
along each spectrum 1 . A datacube representing a single 
SAFARI footprint will therefore be 16x16x1024 pixels in 
size. 

The is the noise commensurate with the sensitivity of SAFARI 
for a lOhr integration time. 
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Fig. 2. A plot of the flux (top row) and redshift (bottom two rows) distributions, for each SED type, of the sources 
that populate the bright-end and burst mode evolution data cubes respectively. The top two rows arc for sources 
with 120 fim flux, Si20m«i > 0.342 mJy and the bottom row is for sources with Si20/im > 1 mJy. Values are for 
a 1 square degree region of sky. Flux is measured as a fraction of the 120 /im confusion limit, <r c = 4.3 mJy 
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A 'truth' catalog is simultaneously generated which 
tracks the location of each source in a datacube, as well 
as the redshift, SED type and 120 /im flux. This allows a 
later comparison between retrieved and input rcdshifts af- 
ter we apply our redshift determination algorithm. This is 
the same datacube and 'truth' catalog generation method 
as used by Clements et al. (2007). 

Datacubes are generated for both burst mode and 
bright-end evolutionary models. For each evolutionary 
model 100 different datacubes are investigated (i.e. taken 
from different regions of the larger 1 square degree cube) 
in order to account for variance in the random spatial 
distributions, each of a size equivalent to 1 SAFARI FoV 
(i.e. a total area of 2'x2'). For each cube the noise along 
the spectra is added by creating 10 differently seeded ran- 
domly generated Gaussian noise arrays in order to account 
for the variance of results due to noise. We therefore have 
a total of 1000 datacubes for each evolutionary model. By 
way of comparison, a figure of a c (X = 120 fim) = 4.3 mJy 
is adopted for the confusion limit, based on Dole et al. 
(2004). 

2.3. Detecting Sources and Extracting Their Redshifts 

Redshifts are determined using a pseudo-cross- 
correlation (PCC) method and stored - however, first the 
presence of a source must be confirmed. This is done by 
checking the 120 ^m flux, S^o^m, of each spatial pixel 
in the cube against a limiting value. This action is per- 
formed before any other operations take place, thereby 
saving processing time on analyzing non-existent or too- 
faint sources. We assume that if any spatial pixel has 
S 120 fj, m > 0.342 mJy (which is equivalent to the la noise 
in a single spectral pixel for 10 hours of integration) then 
a source is present and we attempt to determine its red- 
shift. This value was chosen as our continuum cutoff value 
as empirically it is found that if a source has 120 fxta. flux 
less than this then typically the emission lines are too faint 
to reliably use with our method. A flow diagram illustrat- 
ing the sequence of this method is shown in figure 3. 

Each spectrum is prcprocessed prior to redshift deter- 
mination in the following way: a) Each spectrum is fit 
with a fourth order polynomial which is taken to repre- 
sent the continuum emission of the source; b) the 5i20/im 
value of the spectrum is checked against a limiting value: 
if 5i20mj« < the limiting value, the source is considered 
too faint to determine its redshift, and no further analysis 
is conducted on the spectrum; c) The polynomial contin- 
uum fit is subtracted from the spectrum, leaving an array 
containing only emission lines and noise (sec figure 4) 

The spectral channels in the aforementioned array con- 
taining the 6 highest continuum subtracted flux levels 
(with S\ > 2a) are then initially considered to be emission 
lines and are taken to be our observed lines in our observed 
line array (OLA). Empirically it is found that if the ra- 
tio of the strongest to weakest (continuum subtracted) line 
fluxes present in OLA is less than 1.5 then these lines most 
likely arise from noise and we therefore consider such an 
array to contain no genuine emission lines. Further anal- 
ysis is only conducted on spectra with higher ratio values 



than this. 

OLA is 1024 spectral channels in extent (i.e. identi- 
cal to the spectral extent of our data cube), and contains 
the continuum-subtracted flux of our 6 observed possible 
emission lines at the appropriate spectral channels corre- 
sponding to the wavelengths of those lines, and zero else- 
where. In the case of an emission line stretching across 
multiple pixels, the line is compressed into a single chan- 
nel and there is a built in tolerance in the redshift fitting 
routine to allow for this. The brightest channel present 
in OLA is assumed to hold a genuine emission line. The 
wavelength of this emission line is then compared to the 
set of template wavelengths compiled from a list of the 
strongest emission lines typically seen in galactic spectra 
in the MIR/FIR. This template line array (TLA) is equal 
to one at the appropriate spectral channels correspond- 
ing to the wavelengths of the template emission lines, and 
zero elsewhere. This comparison then provides a table of 
possible rcdshifts the source may lie at. 

If the strongest line in OLA lies at brightest, an d 
our template emission lines lie at wavelengths Xtempiatej , 
where j is the indexing of the line in our template (eg. for 
[OIII]@51.82 fj,m j = 1, for [NIII]@57.32 fj,m j = 2 etc.), 
then an array of possible redshifts at which the source 
may lie can be determined from: 

Xtemplatej 

Where k = 1,2,. . J, where 1 is the total number of lines in 
TLA and thus the total number of possible redshifts. The 
array of template emission lines, TLA, which was used to 
determine the set of possible redshifts is now also used to 
fit to the observed possible emission lines, OLA. Thus 
in order to find the redshift at which TLA most closely 
correlates to OLA, TLA is shifted to each of the possible 
source redshifts, as in; 

Xtemplatej^ — Xtemplatej (1 ~l~ %k) (2) 

We now have I binary rcdshiftcd template line arrays 
(TLAk), each of which are 1024 spectral channels in ex- 
tent and are defined by; 

(TLA)ik(Xi = \template jk ) = 1 (3) 

(TLA) ik (\i^\ templatejk ) = (4) 

where i is the spectral channel corresponding to a given 
wavelength. The strength of the correlation between our 
observed line array and each template line array is given 
by; 

i=1023 

C k = J2 (OLA)i(TLA) ik (5) 

i=Q 

The value of z k which gives the highest value for Ck is 
then assigned as the best estimate of the source redshift. 
In this way the strength of the match depends both on 
the strength of the continuum subtracted emission lines 
that are coincident between the observed spectrum and 
the rcdshiftcd line templates, and the number of matches 
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Fig. 3. A flowchart of the order in which the sources are observed and then our redshift determination method implemented. 
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Fig. 4. Shown in the left panel is a sample spectrum as can be found in our artificial 'skies'. In our 
PCC redshift determination algorithm the continuum of this spectrum is fit with a fourth order poly- 
nomial. This fit is then subtracted from the spectrum, leaving an array containing only emission lines 
and noise. The right panel shows the example spectrum after the continuum fit has been subtracted. 



between lines of the observed and template emission line 
arrays. 

To improve the accuracy of our method, use is made 
of additional, a-priori information, and redshifts which 
produce template/observed array emission line matches 
for strong, commonly observed emission line pairs are 
weighted more heavily. The additional weighting values 
for Cfe when a characteristic pair is found are determined 
empirically to give the most favorable ratio of accurate to 
inaccurate source redshift evaluations 2 3 . Eg. one of the 
strongest line combinations in the MIR is the [SIII]/ [Sill] 
pair at rest wavelengths of 33.42 and 34.82 fim respec- 
tively, both of which are found to be strong in starburst 
galaxies. If a match is made between the observed line 
array and this line pair, then Ck is weighted by an ex- 
tra factor of 2. Two other, weaker, line combinations 
included in our template array are the [OIII]/[NIII] and 
[OIII]/[NII] pairs at rest wavelengths of 51.82, 57.32, 88.36 
and 121.90 (im respectively. If a match is made between 
the observed line array and either of these line pairs, then 
Cfc is weighted by an extra factor of 1.5. 

We include an additional criterion that, as the strengths 
of the members of the [OIII]/[NIII] and [OIII]/[NII] lines 
pairs arc typically comparable, any value of Zk which gives 
a match in the observed line array with one member of 
cither of these pairs, but not its partner, is rejected. This 
relationship holds for all the model spectra used in this 
work, however it may not always hold in practice when 
encountering genuine spectra as observed in extragalactic 

2 An evaluated source redshift is defined as being accurate if 
it differs from the input model redshift by less than 0.1, i.e. 

{^-evaluated ^catalogl < 0.1 

3 Taking as an example three different weighting values: A, B and 
C. A outputs 3 accurate redshifts and inaccurate redshifts. B 
outputs 20 accurate redshifts and 10 inaccurate redshifts. C out- 
puts 15 accurate redshifts and 3 inaccurate redshifts. Of these 
weighting values wc would use C as this outputs a high number 
of accurate redshifts, while limiting the number of inaccurate 
redshifts. 



surveys. 

In contrast to the work described in Clements et al. 
(2007), we can no longer assume that the strongest emis- 
sion line in the spectrum shortward of the continuum emis- 
sion peak (determined from the position of the peak of 
the polynomial fit) is the [OI] line at 63.18 /mi, as we now 
have strong MIR lines present in the spectra. We can, 
however, still assume that the strongest line longward of 
the continuum emission peak is the [CII] line at 157.74 
fxm. This is true as long as the reasonable assumption 
that Tdust > 20 K holds, as [CII] is then typically the only 
strong line longward of the SED peak. Thus, if no single 
redshift is able to map the line template onto the observed 
line array, but the strongest line in the spectrum lies long- 
ward of the SED peak, and has a line to continuum ratio 
> 3 (found empirically to be the lowest value to reliably 
use to identify the [CII] emission line at 157.74 /xm in this 
circumstance) , we assume it to be the [CII] line and, from 
this, calculate a redshift. 

An evaluated redshift is recorded along with the posi- 
tion of the source on the sky. The redshift is referred to 
as the primary redshift and is considered to be that of 
the brightest, or primary, source in the given spatial bin. 
A slightly modified redshift extraction algorithm is then 
run a second-time through the spectrum, to determine 
whether there is a second source of lower flux present; the 
secondary source. When attempting to extract a redshift 
for the secondary source in any spatial bin we first zero 
the lines in OLA which we have already associated with 
the primary source. We then add to OLA the 3 strongest 
lines in the spectrum that have not yet been used in red- 
shift fitting to the spectrum, with flux S\ > 3.5c 4 . This 

4 Eg. defining a line as any spectral channel which has Sx > 2<r, 
if a spectrum contains 20 spectral channels with Sx > 2cr, we 
populate OLA with the brightest 6 of these lines, leaving 14 
unused lines in the spectrum. In order to determine a secondary 
redshift we zero the lines in the OLA which we have already 
associated with the primary source. If, for example, 4 lines in 
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line selection process is used as a) We find empirically that 
6 is the optimum number of emission lines required to ac- 
curately fit a primary redshift (i.e. the maximum amount 
of lines the algorithm is able to use before encountering 
significant degeneracies) b) By definition the secondary 
source is fainter than the primary, thus we expect noise 
to be more of a significant hindrance in redshift fitting 
- therefore we need more stringent requirements on the 
strength of the emission lines used c) 3 is found to be the 
optimum number of lines to add to the previously selected 
lines which are not found to be associated with the pri- 
mary source - if all 6 lines of the initially selected lines 
are found to be associated with the primary source then 3 
emission lines is the minimum number that can be used to 
reliably fit a redshift and any more than this can results in 
significant degeneracies in redshift fitting. The algorithm 
now runs in a manner very similar to before, however (1) 
by definition we expect the continuum-subtracted emis- 
sion lines from the secondary source to be of lower flux 
than those of the primary, and so the requirement on the 
ratio of the strongest to the weakest emission line in the 
observed array is dropped to 1.3, and (2) also by definition 
we are looking at sources of fainter continuum flux, and 
thus are more susceptible to picking up spurious emission 
lines: we therefore no longer weight more heavily redshifts 
which give matches for characteristic MIR/FIR emission 
line pairs, and no longer allow redshifts to be determined 
from the single [CII] emission line. Both these changes in 
the algorithm for secondary as opposed to primary red- 
shift determination are implemented as empirically they 
are found to give the most reliable results. 

The PCC method works by moving through each spec- 
trum in the cube and cross-correlating MIR/FIR line tem- 
plates at a discrete rather than continuous set of redshifts. 
This means that the number of discrete redshifts fitted is 
far less than the number of redshift steps required when 
fitting continuously over a given range. This means we 
are using less of the spectrum and thus are less likely to 
encounter a pixel with a high noise level which could be 
miss-identified as an emission line. We are therefore able 
to use weaker emission lines for redshift determination, 
which enables the algorithm to probe more deeply into 
the noise. This is particularly powerful when looking at 
secondary sources which typically will have weaker line 
fluxes than primary sources 
2.3.1. FIR Emission Lines Only 

As a first test of our automated redshift determination 
we compare how efficiently our method works in compari- 
son to that described in Clements et al. (2007). Clements 
at al. made use of spectra containing FIR emission lines 
only, therefore a direct comparison of this original method, 
and our PCC method as described in section 2.3 can only 

the OLA contributed to the highest value of Cfc then these are 
zeroed, leaving 2 lines remaining in OLA which can be used to 
determine a secondary redshift. Of the 14 remaining unused 
lines above 2cr in the spectrum, 10 of these for example may 
have fluxes S\ > 3.5a, the 3 brightest of which we add to OLA. 
The OLA used to attempt to determine a secondary redshift 
therefore contains a total of 5 lines. 



be made when using exactly the same spectra. To do this 
we use a variant of our PCC method which uses FIR emis- 
sion lines only. The lines are listed in table 2. In addition, 
only the [OIII]/[NIII] and [OIII]/[NII] pairs at 51.82/57.32 
and 88.36/121.90 (im respectively are used to additionally 
weight C'k , as these are the only characteristic pairs which 
lie in the FIR waveband as defined by Clements et al. 
2.3.2. MIR and FIR Emission Lines 

The general version of our method by default uses both 
MIR and FIR lines in our template array, which are listed 
in table 3. The version of the method also allows us to 
make use of all of our characteristic emission line pairs. 

3. Results 

All results given in this section were determined by tak- 
ing the averaged values for each of 100 FoVs (x 10 differ- 
ently seeded noise arrays) for both bright-end and burst 
mode evolution. When quoting results from this work we 
give the number of accurate redshifts retrieved as a per- 
centage of the total number of sources within a given flux 
and redshift range in the datacube. However the number 
of inaccurate redshifts are given as a percentage of the 
total number of redshifts output by the algorithm under 
the same constraints 5 . 

3.1. Results From Method of Clements et al. (2007) 

The work described in Clements et al. (2007) made use 
of spectra containing FIR emission lines only, and results 
were only determined for the burst mode evolution model. 
They found that all sources with redshifts at z < 2.5 with 
Si20fim > 1 mjy could be retrieved. Redshifts higher than 
this could not be determined, as beyond z = 2.5 the [OI] 
and [CII] lines are redshiftcd out of the SAFARI pass- 
band. Lower flux sources, as much as ^10 times fainter 
than the traditional continuum confusion limit were also 
retrieved, albeit with lower efficiency. 

3.2. Analysis Using FIR Emission Lines Only 

3.2.1. Primary Sources 

Using the FIR emission line only version of the method 
outlined in section 2.3.1 we find that we recover accu- 
rate redshifts for 85(46)% of all primary sources with 
Si20fmi — 1 rnjy, with 5(10)% of all primary redshifts out- 
put by the algorithm under the same constraints being 
inaccurate. The aforementioned recovery values are given 
as a percentage of all sources, including those at z > 2.5, 
however by definition this method is unable to retrieve 
redshifts for sources at z > 2.5. Beyond z = 2.5 the redshift 
recovery rate drops to zero for both evolutionary models. 
For sources with 5i20^m > 1 rnJy lying at z < 2.5 the recov- 
ery of redshifts for primary sources is ~85(64)%. Under 
the same constraints Clements et al. (2007) retrieved ac- 

5 Eg. a datacube has 100 sources within the specified redshift and 
flux range and the algorithm outputs 50 redshifts; 40 accurate 
and 10 inaccurate. The results statement would then be given 
in the following form; we retrieve 40% of sources in our flux and 
redshift range accurately, with 20% of redshifts output by the 
algorithm under the same constraints being inaccurate. 
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Emission Line OIII NIII 01 OIII Nil 01 CII 
Wavelength (pm) 51.82 57.32 63.18 88.36 121.90 145.53 157.74 

Table 2. Lines used in source redshift determination through template fitting for the FIR emission line only method. 



Emission Line Nell SIII SIII Sill OIII NIII 01 OIII Nil 01 CII 
Wavelength (/mi) 12.81 18.71 33.42 34.82 51.82 57.32 63.18 88.36 121.90 145.53 157.74 

Table 3. Lines used in source redshift determination through template fitting for the method using both FIR and MIR. emission lines. 



curate redshifts for 100% of the sources. However in that 
work redshifts were determined manually for each source 
whereas in this work they are determined automatically. 
Unlike Clements ct al. we are unable to assess each spec- 
trum on a case by case basis. Therefore in attempting to 
minimize the number of incorrect redshifts output by the 
algorithm we are forced to limit the maximum possible ef- 
ficiency of the accurate redshift recovery of the algorithm. 

By dropping our 120 /tm flux cutoff to S\20^m > 
0.342 mJy, we find that we recover accurate redshifts for 
36(33)% of all primary sources with Si20/j.m > 0.342 mJy, 
with 12(13)% of all redshifts output by the algorithm un- 
der the same constraints being inaccurate. Taking into 
account only sources at z < 2.5 our accurate recovery be- 
comes 45(39)%. Despite this being a lower recovery per- 
centage than for Si20nm > 1 m Jy the majority of sources 
have lower fluxes than this. As such by dropping the cut- 
off flux from 1 mJy to 0.342 mJy we have increased the 
total number of accurate redshifts recovered by 22(15)%. 
3.2.2. Secondary Sources 

Using our FIR emission line only PCC method we re- 
trieve accurate redshifts for ^30(27)% of all secondary 
sources when using a flux cutoff of S^o^m > 0.342 mJy 
with ~6(4)% of all secondary redshifts output by the al- 
gorithm under the same constraints being inaccurate. 

3.3. Extending Analysis to MIR 

By default, the full version of our PCC redshift de- 
termination method uses both MIR and FIR emission 
lines, and is that version which would be applied to real 
data. It is now possible to investigate sources at redshifts 
z > 2.5 which account for 5(26)% of the population with 
Si20iim > 0.342 mJy, of the evolutionary models. Accurate 
redshifts are retrieved for ^75% of primary sources with 
S '120 fj,m > 1 rnjy, for both bright-end and burst mode evo- 
lution, with respectively 6(8)% of all primary redshifts 
output by the algorithm under the same constraints be- 
ing inaccurate. We see here that using both MIR and 
FIR emission lines is less efficient in recovering accurate 
redshifts than using FIR emission lines only at redshifts 
z < 2.5. This is due to the fact the ability to resolve red- 
shifts z > 2.5 is not as advantageous here as most sources 
with S\20ftm > 1 m Jy lie at redshifts z < 2.5 in both evo- 
lution models, however we are still subject to the disad- 
vantage of significantly more degeneracies in redshift fit- 
ting due to a higher number of emission lines. Dropping 
the 120 /urn flux cutoff, 38(54)% of all primary sources 
with Si20fim > 0.342 mJy arc retrieved (see figure 5), with 
14(9)% of all redshifts output by the algorithm under the 



same constraints being inaccurate. 
3.3.1. Secondary Sources 

Employing a 120 /im cutoff of 0.342 mJy, wc find that we 
are able to determine accurate redshifts for ~38(29)% of 
secondary sources with Si20fim > 0.342 mJy (see figure 6), 
with 11(18)% of all redshifts output by the algorithm un- 
der the same constraints being inaccurate. 

4. Discussion 

Wc find that using the PCC method as described in 
this work, we can recover redshifts for sources as much as 
10 times below the traditional continuum confusion limit. 
We also find we are able to successfully retrieve redshifts 
from line confused spectra caused by multiple sources sep- 
arated by smaller scales than the size of the beam on the 
sky. However, using the burst mode and bright-end evolu- 
tionary models of Pearson (2005), Pearson et al. (2007) 
and Pearson & Khan (2009), we find no spatial bins con- 
taining more than two sources with Si20ftm > 0.342 mJy 
and therefore we do not test whether or not our method 
would be able to disentangle more sources than this. 

Shown in figure 7 is the cumulative recovery (fraction 
of sources accurately recovered with S'120/im less than that 
defined by the x-axis) efficiency with increasing flux. All 
source populations (the primary and secondary sources of 
both bright-end and burst mode evolution) have a strongly 
increasing cumulative recovery fraction at low fluxes which 
then levels off at fluxes higher than the confusion limit. 
This is because the bulk of their populations lie at fluxes 
fainter than the confusion limit. 

The bright-end evolution model has a larger fraction of 
low flux sources, whereas the burst mode evolution model 
has a similar population of low flux sources, but also a 
larger number of brighter sources. The efficiency of our 
method falls off with decreasing flux (as illustrated in fig- 
ure 5), therefore we retrieve the redshifts for the burst 
mode evolution model more efficiently than for the bright- 
end evolution model. 

At higher fluxes the redshift determination efficiency 
for all populations goes to 100%, however at these higher 
fluxes there are relatively fewer sources. Even at 10% 
of the confusion limit, we are still retrieving redshifts for 
~10% of sources, and given that the number of sources at 
these fluxes is so large we are gaining information about a 
significantly larger number of galaxies than would be the 
case if we were confusion limited. 

The percentage of retrieved redshifts which are in er- 
ror is higher for secondary sources than for the primary 




Fig. 5. Plots of the 120 fim continuum flux distribution (with flux measured as a fraction of the 120 fim con- 
tinuum confusion limit, cr c (A = 120^m) = 4.3 mjy) for both the total input primary sources (dashed) and pri- 
mary sources with accurately determined redshifts (solid). Shown in the left and right hand panels are the re- 
sults for the bright-end and burst mode evolution models respectively. Results arc for 100 SAFARI FoVs. 
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Fig. 6. Plots of the 120 fim continuum flux distribution (with flux measured as a fraction of the 120 fim con- 
tinuum confusion limit, <r c (A = I20fim) = 4.3 mjy) for both the total input secondary sources (dashed) and sec- 
ondary sources with accurately determined redshifts (solid). Shown in the left and right hand panels are the 
results for the bright-end and burst mode evolution models respectively. Results are for 100 SAFARI FoVs. 



sources. This is to be expected because by definition the 
secondary sources are fainter than the primary. We are 
thus more likely to confuse noise with genuine emission 
lines. Another difficulty encountered when attempting to 
retrieve secondary redshifts is that we are unable to set 
a continuum flux limit to consider the source to be vi- 
able. Therefore we are most often attempting to retrieve 
redshifts from either very low flux or non-existent sec- 
ondary galaxies, significantly increasing the likelihood of 
retrieving an erroneous redshift. This problem is lessened 
somewhat by the introduction of more stringent criteria 
in other areas when attempting to retrieve secondary red- 
shifts, but it remains the cause of a large percentage of 
the inaccurate secondary redshift recoveries. 

We find that by employing a cross correlation method 
for redshift determination, where we are only consider- 
ing a very discrete set of redshift possibilities, we are able 



to reduce the Si20y.m cutoff value without drastically de- 
creasing the efficiency of the method and thus are able 
to dig deeper into the instrumental noise. Employing the 
full MIR and FIR emission lines version of our method 
we find that we can retrieve accurate redshifts for a total 
of 38(52)% of all (i.e. both the primary and secondary 
sources) sources with Si20/xm > 0.342 mjy for bright-end 
and burst mode evolution respectively. Their respective 
frequency of occurrence of inaccurate redshifts as a per- 
centage of all redshifts output by the algorithm under the 
same constraints is 15(10)%. 

A limitation of this method is the possibility of sources 
with anomalous line strengths which can, for example, re- 
sult in a miss-identification of the [CII] line. The possibil- 
ity of such occurrences has not been taken in to account 
here. Further improvement to our PCC method would 
aim to decrease the total number of erroneous redshifts 



12 




Fig. 7. Plots of the cumulative fractional recovery of sources (fraction of sources accurately recovered with Si20/j.m less than 
that defined by the x-axis) with increasing 120 (im flux, measured as a fraction of the traditional 120 /jm continuum confusion 
limit, <r c (A = 120/j.m) = 4.3 mjy. The left panel shows our results for primary sources and the right panel shows our results for 
secondary sources. The results for burst mode evolution are plotted as a solid line, and for bright-end evolution as a dashed line. 



output for reasons such as this. 

Additional spectral features that can be used for finding 
redshift are those associated with the PAHs (Polycyclic 
Aromatic Hydrocarbons). These features are much 
broader than the emission lines we have been using, thus 
observations could be made at lower spectral resolution, 
with a corresponding increase in instrument sensitivity. 
In this work we have modeled angular resolution as being 
constant over the full SAFARI band. If we were to em- 
ploy a more realistic model for angular resolution we may 
observe sources which arc clustered in a single spatial bin 
in our current 8" binning, as being separable (visible in 
different spatial bins) at shorter wavelengths. This could 
significantly decrease the fraction of erroneous redshifts 
output for secondary sources. 

As discussed earlier, the sources from the bright-end 
evolution model are mostly at low fluxes whereas the burst 
mode has in addition a smaller population of high flux 
sources. This accounts for the difference in our recovery 
rates for the two models. It is therefore important to 
further test the efficacy of our technique, and of this deep 
observing mode of SAFARI spectral imaging on a number 
of different evolutionary models. 

Using a Kolmogorov-Smirnov test we are able to deter- 
mine whether a set of recovered redshifts are from differ- 
ent parent populations. Thus by comparing redshift dis- 
tributions recovered from datacubes of different sizes we 
can determine what area of sky is required for SAFARI 
to be able to reliably distinguish between different evo- 
lutionary models. We compare the maximum deviation 
between the cumulative distributions (Wall 1996) of the 
redshifts output from our PCC method for skies popu- 
lated with the burst mode evolution model and with the 
bright-end evolution model. We perform this comparison 
for increasing numbers of SPICA FoVs. Using our PCC 
redshift determination method we find we are able to reli- 
ably distinguish (probability that the sources from the two 
models are drawn from the same distribution, P =0.01%) 



between the bright-end and burst mode evolution models 
with a sky survey area equal to 8 SAFARI FoVs, each of 
10 hrs integration time. 

5. Conclusions 

We have found that our PCC redshift determination 
method is capable of resolving sources (i.e. determining a 
unique redshift) more than an order of magnitude fainter 
than the traditional continuum confusion limit, however 
the efficacy of our method is higher for brighter sources. 
In this work we have used the PCC method on models 
based around the SAFARI instrument for SPICA, how- 
ever the same technique could be used on any sensitive 
imaging spectrometer. The bright-end and burst mode 
evolution models include sources up to redshifts of z ~4 
and 5 respectively. At these redshifts we are still able to 
determine redshifts for sources using the PCC method. 
We have not yet tested to see at which redshift the PCC 
method begins to fail, this may be the subject of future 
work. 

The evolutionary models we have investigated in this 
work have the bulk of their populations at fluxes fainter 
than the confusion limit. By employing our PCC method 
we are therefore greatly increasing the number of galaxies 
that we are able to uniquely identify in any single obser- 
vation. This presents us with a better statistical sample 
with which to compare observed source counts and red- 
shift distributions with those presented in different evo- 
lutionary models. We therefore have the potential to re- 
liably distinguish different evolutionary models and our 
observations with much smaller area surveys, and there- 
fore within shorter observing times. 

Future work should include the following: 1) Use of 
the PAH features to identify sources and determine their 
redshifts. This will allow us to decrease our spectral 
resolution, thus increasing instrument sensitivity. 2) 
Investigation of the viability of taking into account sources 
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with atypical line strengths. 3) More realistic angular res- 
olution modeling where spatial resolution varies across the 
waveband. 4) Investigation of the efficiency of the method 
when implemented on a wider range of evolutionary mod- 
els. 5) A more quantitative analysis of where our ability 
to retrieve redshifts from single and combined spectra be- 
gins to break down is being conducted. 6) It should also 
be noted that since the work described in this paper was 
conducted the technical specifications of SAFARI have 
changed somewhat (eg. waveband, sensitivity), therefore 
the results should be re-checked with more up to date 
modeling of the SAFARI instrument. These topics are 
currently being investigated. 
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